home *** CD-ROM | disk | FTP | other *** search
/ Tech Arsenal 1 / Tech Arsenal (Arsenal Computer).ISO / tek-01 / ddj9304.zip / WAVELET.ZIP / WAVE_BLD.C < prev    next >
Text File  |  1992-05-18  |  9KB  |  331 lines
  1. /* WAVE_BLD.C */
  2.  
  3. #define WAVE_BLD
  4.  
  5. #include <math.h>
  6. #include <alloc.h>
  7. #include "wave_bld.h"
  8.  
  9. typedef enum {DECOMP, RECON} wavetype;
  10.  
  11. double **BuildWaveStorage(int inlength, int *levels);
  12. void DestroyWaveStorage(double **tree, int levels);
  13. void WaveletCoeffs(double alpha, double beta, double *wavecoeffs);
  14. char MakeWaveletFilters(double *wavecoeffs, double *Lfilter,
  15.                                                             double *Hfilter, wavetype transform);
  16. double DotpEven(double *data, double *filter, char filtlength);
  17. double DotpOdd(double *data, double *filter, char filtlength);
  18. void ConvolveInt2(double *input_sequence, int inp_length, double *filter,
  19.                 char filtlen, char sum_output, double *output_sequence);
  20. void ScalingFuncRegeneration(double **InData, int Inlength,
  21.                 double *Lfilter, char filtlen, char levels, double *OutData);
  22. void WaveletRegeneration(double **InData, int Inlength, double *Lfilter,
  23.                                 double *Hfilter, char filtlen, char levels, double *OutData);
  24. double DataSetMagnitude(double *data, int length);
  25.  
  26. double wavedata[901], wave_coeff[6];
  27. double g[6], gtilda[6], h[6], htilda[6];
  28. double **ExCoeffPhi, **ExCoeffPsi;
  29. int wave_levels;
  30.  
  31. double **BuildWaveStorage(int inlength, int *levels)
  32. {
  33.     double **tree;
  34.     int i, j, lvls;
  35.  
  36.     lvls = *levels;
  37.     /* create decomposition tree */
  38.     tree = (double **) calloc(lvls, sizeof(double *));
  39.     j = inlength;
  40.     for (i = 0; i < lvls; i++)
  41.     {
  42.         j /= 2;
  43.         if (j == 0) /* too many levels requested for available data */
  44.         {           /*  number of transform levels now set to value of i */
  45.             *levels = i;
  46.             continue;
  47.         }
  48.  
  49.         tree[i] = (double *) calloc((j + 5), sizeof(double));
  50.     }
  51.  
  52.     return tree;
  53. }
  54.  
  55.  
  56. void DestroyWaveStorage(double **tree, int levels)
  57. {
  58.     char i;
  59.  
  60.     for (i = levels - 1; i >= 0; i--)
  61.         free(tree[i]);
  62.  
  63.     free(tree);
  64. }
  65.  
  66.  
  67. void WaveletCoeffs(double alpha, double beta, double *wavecoeffs)
  68. {
  69.     double tcosa, tcosb, tsina, tsinb;
  70.     char i;
  71.     /* precalculate cosine of alpha and sine of beta to reduce */
  72.     /* processing time */
  73.     tcosa = cos(alpha);
  74.     tcosb = cos(beta);
  75.     tsina = sin(alpha);
  76.     tsinb = sin(beta);
  77.  
  78.     /* calculate first two wavelet coefficients  a = a(-2) and b = a(-1) */
  79.     wavecoeffs[0] = ((1.0 + tcosa + tsina) * (1.0 - tcosb - tsinb)
  80.                                                                             + 2.0 * tsinb * tcosa) / 4.0;
  81.     wavecoeffs[1] = ((1.0 - tcosa + tsina) * (1.0 + tcosb - tsinb)
  82.                                                                             - 2.0 * tsinb * tcosa) / 4.0;
  83.  
  84.     tcosa = cos(alpha - beta); /* precalculate cosine and sine of alpha */
  85.     tsina = sin(alpha - beta); /* minus beta to reduce processing time */
  86.  
  87.     /* calculate last four wavelet coefficients  c = a(0), d = a(1), */
  88.     /* e = a(2), and f = a(3) */
  89.     wavecoeffs[2] = (1.0 + tcosa + tsina) / 2.0;
  90.     wavecoeffs[3] = (1.0 + tcosa - tsina) / 2.0;
  91.     wavecoeffs[4] = 1 - wavecoeffs[0] - wavecoeffs[2];
  92.     wavecoeffs[5] = 1 - wavecoeffs[1] - wavecoeffs[3];
  93.  
  94.     /* zero out very small coefficient values caused by truncation error */
  95.     for (i = 0; i < 6; i++)
  96.     {
  97.         if (fabs(wavecoeffs[i]) < 1.0e-15)
  98.             wavecoeffs[i] = 0.0;
  99.     }
  100. }
  101.  
  102.  
  103. char MakeWaveletFilters(double *wavecoeffs, double *Lfilter,
  104.                                                             double *Hfilter, wavetype transform)
  105. {
  106.     char i, j, k, filterlength;
  107.  
  108.     i = 0; /* find the first non-zero wavelet coefficient */
  109.     while(wavecoeffs[i] == 0.0)
  110.         i++;
  111.  
  112.     j = 5; /* find the last non-zero wavelet coefficient */
  113.     while(wavecoeffs[j] == 0.0)
  114.         j--;
  115.  
  116.     /* form the decomposition filters h~ and g~ or the reconstruction
  117.          filters h and g.  the division by 2 in the construction
  118.          of the decomposition filters is for normalization */
  119.     filterlength = j - i + 1;
  120.     for(k = 0; k < filterlength; k++)
  121.     {
  122.         if (transform == DECOMP)
  123.         {
  124.             Lfilter[k] = wavecoeffs[j--] / 2.0;
  125.             Hfilter[k] = (double) (((i++ & 0x01) * 2) - 1) * wavecoeffs[i] / 2.0;
  126.         }
  127.         else
  128.         {
  129.             Lfilter[k] = wavecoeffs[i++];
  130.             Hfilter[k] = (double) (((j-- & 0x01) * 2) - 1) * wavecoeffs[j];
  131.         }
  132.     }
  133.  
  134.     while (k < 6) /* clear out the additional array locations, if any */
  135.     {
  136.         Lfilter[k] = 0.0;
  137.         Hfilter[k++] = 0.0;
  138.     }
  139.  
  140.     return filterlength;
  141. }
  142.  
  143.  
  144. double DotpEven(double *data, double *filter, char filtlen)
  145. {
  146.     char i;
  147.     double sum;
  148.  
  149.     sum = 0.0;
  150.     for (i = 0; i < filtlen; i += 2)
  151.         sum += *data-- * filter[i]; /* decreasing data pointer is moving */
  152.                                                                 /* backward in time */
  153.     return sum;
  154. }
  155.  
  156.  
  157. double DotpOdd(double *data, double *filter, char filtlen)
  158. {
  159.     char i;
  160.     double sum;
  161.  
  162.     sum = 0.0;
  163.     for (i = 1; i < filtlen; i += 2)
  164.         sum += *data-- * filter[i];    /* decreasing data pointer is moving */
  165.                                                                 /* backward in time */
  166.     return sum;
  167. }
  168.  
  169.  
  170. void ConvolveInt2(double *input_sequence, int inp_length, double *filter,
  171.                     char filtlen, char sum_output, double *output_sequence)
  172.     /* insert zeros between each element of the input sequence and
  173.        convolve with the filter to interpolate the data */
  174. {
  175.     int i;
  176.  
  177.     if (sum_output) /* summation with previous convolution if true */
  178.     {
  179.             /* every other dot product interpolates the data */
  180.         for(i = (filtlen / 2) - 1; i < inp_length + filtlen - 2; i++)
  181.         {
  182.             *output_sequence++ += DotpOdd(input_sequence + i, filter, filtlen);
  183.             *output_sequence++ += DotpEven(input_sequence + i + 1, filter, filtlen);
  184.         }
  185.  
  186.         *output_sequence++ += DotpOdd(input_sequence + i, filter, filtlen);
  187.     }
  188.     else /* first convolution of pair if false */
  189.     {
  190.             /* every other dot product interpolates the data */
  191.         for(i = (filtlen / 2) - 1; i < inp_length + filtlen - 2; i++)
  192.         {
  193.             *output_sequence++ = DotpOdd(input_sequence + i, filter, filtlen);
  194.             *output_sequence++ = DotpEven(input_sequence + i + 1, filter, filtlen);
  195.         }
  196.  
  197.         *output_sequence++ = DotpOdd(input_sequence + i, filter, filtlen);
  198.     }
  199. }
  200.  
  201.  
  202. void ScalingFuncRegeneration(double **InData, int Inlength,
  203.                 double *Lfilter, char filtlen, char levels, double *OutData)
  204.     /* assumes that the input data has (2^levels) data points per unit
  205.        interval */
  206. {
  207.     double *Output;
  208.     char i;
  209.  
  210.     Inlength = Inlength >> levels;
  211.     /* destination of all but last branch reconstruction is the next
  212.          higher intermediate approximation */
  213.     for (i = levels - 1; i > 0; i--)
  214.     {
  215.         Output = InData[i - 1];
  216.         ConvolveInt2(InData[i], Inlength, Lfilter, filtlen, 0, Output);
  217.         Inlength *= 2;
  218.     }
  219.  
  220.     /* destination of the last branch reconstruction is the output data */
  221.     ConvolveInt2(InData[0], Inlength, Lfilter, filtlen, 0, OutData);
  222. }
  223.  
  224.  
  225. void WaveletRegeneration(double **InData, int Inlength, double *Lfilter,
  226.                                 double *Hfilter, char filtlen, char levels, double *OutData)
  227.     /* assumes the input data has (2^levels) data points per unit interval */
  228. {
  229.     double *Output;
  230.     char i;
  231.  
  232.     Inlength = Inlength >> levels;
  233.     ConvolveInt2(InData[levels - 1], Inlength, Hfilter, filtlen,
  234.                                                                                             0, InData[levels - 2]);
  235.     Inlength *= 2;
  236.     /* destination of all but last branch reconstruction is the next
  237.          higher intermediate approximation */
  238.     for (i = levels - 2; i > 0; i--)
  239.     {
  240.         Output = InData[i - 1];
  241.         ConvolveInt2(InData[i], Inlength, Lfilter, filtlen, 0, Output);
  242.         Inlength *= 2;
  243.     }
  244.  
  245.     /* destination of the last branch reconstruction is the output data */
  246.     ConvolveInt2(InData[0], Inlength, Lfilter, filtlen, 0, OutData);
  247. }
  248.  
  249.  
  250. double DataSetMagnitude(double *data, int length)
  251. {
  252.     int i;
  253.     double newval, maxval;
  254.  
  255.     if (data[0] < 0) /* initial seed for largest magnitude of data set */
  256.         maxval = -data[0];
  257.     else
  258.         maxval = data[0];
  259.  
  260.     for (i = 1; i < length; i++)
  261.     {                /* determine largest magnitude of data set */
  262.         if (data[i] < 0)
  263.             newval = -data[i];
  264.         else
  265.             newval = data[i];
  266.  
  267.         if (newval > maxval)
  268.             maxval = newval;
  269.     }
  270.  
  271.     if (maxval == 0.0) /* default value if all data points are zero */
  272.         maxval = 1.0;
  273.  
  274.     return maxval;
  275. }
  276.  
  277.  
  278. void SetupWaveletConstruction(void)
  279. {
  280.     wave_levels = 6;
  281.     ExCoeffPhi = BuildWaveStorage(448, &wave_levels);
  282.     ExCoeffPsi = BuildWaveStorage(448, &wave_levels);
  283.     ExCoeffPhi[5][5] = 1.0;   /* reconstruction from a single impulse ... */
  284.     ExCoeffPsi[5][5] = 1.0;   /* yields the scaling or wavelet function */
  285.     PhiData = wavedata;       /* both functions are in one array */
  286.     PsiData = &wavedata[448];
  287. }
  288.  
  289.  
  290. void BuildDSPfilterArray(double alpha, double beta, float *DSPfilters)
  291. {
  292.     int i, j, filtlen;
  293.  
  294.     /* determine wavelet coefficients for DSP */
  295.     WaveletCoeffs(alpha, beta, wave_coeff);
  296.     MakeWaveletFilters(wave_coeff, htilda, gtilda, DECOMP);
  297.     /* place filters in array for download to DSP board */
  298.     i = 0;
  299.     for (j = 0; j < 6; j++)
  300.         DSPfilters[i++] = (float) gtilda[j];
  301.  
  302.     for (j = 0; j < 6; j++)
  303.         DSPfilters[i++] = (float) htilda[j];
  304.  
  305.     /* determine wavelet coefficients for function calculations */
  306.     filtlen = MakeWaveletFilters(wave_coeff, h, g, RECON);
  307.  
  308.     for (i = 896; i < 901; i++) /* zeroo out end effects */
  309.         wavedata[i] = 0.0;
  310.  
  311.     ScalingFuncRegeneration(ExCoeffPhi, 448, h, filtlen, 6, PhiData);
  312.     WaveletRegeneration(ExCoeffPsi, 448, h, g, filtlen, 6, PsiData);
  313.     magnitude = DataSetMagnitude(wavedata, 901);
  314.  
  315.     /* depending upon length of the filters,
  316.          the functions will be offset by a fixed amount */
  317.     switch(filtlen)
  318.     {
  319.         case 2:
  320.             offset = 128;
  321.             break;
  322.  
  323.         case 4:
  324.             offset = 64;
  325.             break;
  326.  
  327.         case 6:
  328.             offset = 0;
  329.     }
  330. }
  331.